P'iWk:  reporting  burdon  for  ; 

n  liinTalning  the  data  needec  1 1 

tor  reducing  !h)«  burdarr,  lo 
!ha  Office  of  Management  i 

1 .  Agency  Use  Only  (Leave  blank).  j  2.  Report  Dale. 

July  1992 


4.  Title  and  Subtitle. 

Least-squares  time-delay  estimation  for  transient  signals  in 


6.  Auihor(s). 

R,  J,  Vaccaro*.  C.  S.  Ramalingam*,  D.  W.  Tufts',  and  R.  L. 


7.  Performing  Organization  Name^s)  and  Address(es). 
Naval  Research  Laboratory 
Acoustics  Division 

Stennis  Space  Center,  MS  39529-5004 


r  Approved 
OBMNa.  ‘y/04  0W8 


kfding  the  linm  tor  reviewing  intiru'  tf'jns.  so.trc'i,<ig  e»'5fi*  9  '.TfS  g  ’> 

jardinglhisdufdenof  any  ,T’;urdM)eciofth!«ccil9cfiofiof  .thon.sr'ck/d  '^s 

kJ  Reports.  121s  Je(tu;son  Davis  Highway,  Suse  Artu-T’on,  VA  4 302  ; 

ix _ _ _  — 

3.  Report  Type  and  Dales  Covered.  ^  -  -j 

Final  -  Journal  Article  \ 


S.  Funding  Numbers.  V., ^ 

a  multipath  environment  Program  Element  No  0601153N 


Project  No. 

Task  No 
Accession  No. 
Work  Unit  t.’ 


OHO 

UN257015 
:  ■^4doC 


B.  Performlnn  -''’nizatlcn 
Report  Nuni*’t'- 

JA  244:048:91 


9.  Sponsoring/Monitoring  Agency  Name(s)  and  Address(es), 
Naval  Research  Laboratory 

Center  for  Environmental  Acoustics  p 

Stennis  Space  Center.  MS  39529-5004 


ELECTt 
MAY  17  1993 


11.  Supplementary  Notes. 

Published  in  J.  Acoust.  Soc.  Am. 

'The  University  of  Rhode  Island.  Kingston,  Rl  39529-0288 


10.  Sponsorlng/Moi.'iorlng  Agency 
Report  Number. 

JA  244:048:91 


12a.  Distrlbutlon/Availablllty  Statement. 

Approved  for  public  release;  distribution  is  unlimited. 


12b.  Distribution  Cede, 


13.  Abstract  (Maximum  200  words). 

Theproblemof  estimating  the  arrival  times  of  overlapping  ocean-acouslic  sign  'is  from  a  noisy  received  waveh'rmlh  ii 
of  scaled  and  delayed  replicas  of  a  deterministic  transient  signal  is  considered  I'  is  assunu-d  that  the  transmit!’  i  sic  .ai  t  '  n, 
number  of  paths  in  the  muHipath  environment  are  known,  and  an  algorithm  is  d  /eloped  Hiat  gives  least  squ-s’es  pstima'e?  c* 
the  amplitude  and  time  delay  of  each  path.  A  method  is  given  to  ensure  that  the  riiobal  minimum  of  the  error  surface  is  found  in 
spite  of  the  existence  of  numerous  local  minima.  The  algorithm  is  then  extended  ’n  the  case  in  which  the  transm-tted  signal  is  not 
known  precisely,  but  is  assumed  to  belong  to  a  parametric  class  of  signals.  Tfie  extended  algorithm  additionally  obtains  the 
parameters  that  characterize  the  transmitted  signal.  The  algorithm  is  demonstrated  on  the  class  of  signals  consisting  of  gated 
sinusoids,  using  both  simulated  and  experimental  data. 


93 


0 


93-10903 

lilillllllllll'IlvMn 


14.  Sub)«ct  Term*. 

Transients,  distributed  sensors,  coherence,  detection,  classification 


17.  Security  Clacsificatlon 
of  Report. 

Unclassified 


NSN  7540  01  280  5500 


18,  Security  Classification 
of  This  Page. 

Unclassified 


19.  Security  rtlasslflcatlon 
of  Absirnci, 

Unclassilied 


15.  Number  of  Pages. 

9 


16.  Price  Code. 


[  20.  Limitation  nf  Abstract. 

I  SAR 


'"■  (  orm  ?98  fRov  ?  0?) 

r  .W.nbF>«i  >  V  ANSI  Sid  fS 
10? 


Least-squares  time-delay  estimation  for  transient  signals 
in  a  multipath  environment 

R.J.  Vaccaro,  C.  S.  Ramalingam.andD.  W  Tufts 

Department  of  Electrical  Engineering,  The  University  of  Rhode  hland,  Kingston,  Rhode  Island  02SSI 

I  R.  L  Field 
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I  (Received  ISJune  1991;revised  1  October  1991; accepted 9 February  1992) 

The  problem  of  estimating  the  arrival  times  of  overlapping  ocean-acoustic  signals  from  a  noisy 
received  waveform  that  consists  of  scaled  and  delayed  replicas  of  a  deterministic  transient 
,  signal  is  considered.  It  is  assumed  that  the  transmitted  signal  and  the  number  of  paths  in  the 

multipath  environment  are  known,  and  an  algorithm  is  developed  that  gives  least-squares 
i  estimates  of  the  amplitude  and  time  delay  of  each  path.  A  method  is  given  to  ensure  that  the 

global  minimum  of  the  error  surface  is  found  in  spite  of  the  existence  of  numerous  local 
minima.  I'he  algorithm  is  then  extended  to  the  case  in  which  the  transmitted  signal  is  not 
known  precisely,  but  is  assumed  to  belong  to  a  parametric  class  of  signals.  The  extended 
algorithm  additionally  obtains  the  parameters  that  chrtiacterlzc  the  Uansiniited  signal.  1  he 
algorithm  is  demonstrated  on  the  class  of  signals  consisting  of  gated  sinusoids,  using  both 
simulated  and  experimental  data. 

PACS  numbers:  43.60.Gk,  43.30Wi,  43.30.Ma 


INTRODUCTION 

Time  delay  estimation  is  a  well-known  problem  occur¬ 
ring  frequently  in  the  fields  of  sonar,  radar,  and  geophysics. 
In  this  problem,  the  received  waveform  consists  of  delayed 
and  scaled  replicas  of  the  transmitted  signal.  This  is  the  re¬ 
sult  of  multiple  reflections  and  attenuation  of  the  signal  in 
the  channel. 

The  received  waveform  r{t)  can  be  described  math¬ 
ematically  as 

SI 

r{t)~  ^  ai,s(t  -  T^)  +  nU),  0<r<r,  (1) 

where  s(t)  is  the  transmitted  signal,  is  the  attenuation 
factor  for  path  k,  is  the  time  delay  for  path  k,  M  is  the 
number  of  different  paths,  and  n(r)  is  the  inevitable  noise 
component  corrupting  the  received  signal.  We  say  that  s(t) 
is  a  deterministic  signal,  which  means  that  we  do  not  assume 
any  statistical  properties  for  the  signal  such  as  its  power 
spectra!  density.  The  signal  waveform  itself  is  used  by  the 
algorithms  developed  in  this  paper.  For  the  purposes  of  this 
paper,  we  define  a  transient  signal  to  be  a  signal  whose  dura¬ 
tion  is  less  than  2  s. 

Even  though  Eq.  ( 1 )  has  been  stated  in  continuous  time, 
any  practical  implementation  using  digital  processing  tech¬ 
niques  would  deal  only  with  discrete-time  samples.  Hence, 
in  what  follows  we  will  consider  only  discrete-time  signals. 
Further,  we  will  assume  that  the  noise  is  white  Gau.ssian.  We 
will  also  assume  that  the  number  of  paths  A/ is  known.  This 
latter  assumption  is  not  as  restrictive  as  it  appears  to  be,  since 
in  many  cases  the  number  of  paths  can  be  determined  quite 
reasonably  from  the  bathymetry  of  the  channel,  the  approxi¬ 
mate  locations  of  the  transmitter  and  receiver,  the  sound 
velocity  profile,  and  a  propagation  model.' 


The  classical  method  for  estimating  the  times  of  arrival 
is  correlating  the  received  waveform  with  the  transmitted 
waveform.^  The  peaks  in  the  correlator  output  give  the  esti¬ 
mates  of  the  arrival  times.  It  can  be  shown  that  if  the  signals 
are  separated  in  time  by  more  than  the  duration  of  the  signal 
autocorrelation  function,  the  correlator  is  equivalent  to  the 
maximum-likelihood  estimator.^  Other  approaches  are  giv¬ 
en  in  Refs.  4,  5,  6,  and  7. 

A  completely  different  approach  was  recently  proposed 
by  Kirsteins.* ''  The  basic  idea  of  this  approach  is  to  view  *he 
problem  in  the  frequency  domain.  Since  a  delay  in  ihe  time 
domain  is  equivalent  to  multiplication  by  an  exponential  in 
the  frequency  domain,  the  corresponding  frequency  domain 
problem  is  one  of  fitting  weighted  complex  exponentials  to 
the  spectrum  of  the  received  signal.  Utilizing  an  iterative 
method  of  fitting  complex  exponentials  as  in  Refs.  10  and  11, 
this  approach  provides  a  way  of  estimating  the  times  of  arriv¬ 
al.  However,  this  algorithm  does  not  work  if  the  Fourier 
transform  of  the  source  signal  has  any  nulls  (zeros)  in  the 
regions  of  its  spectral  support.  Even  if  a  spectral  sample  does 
not  fall  exactly  on  a  null,  a  small  sample  value  could  cause 
serious  numerical  ill-conditioning  of  the  algorithm.  Hence 
we  cannot  apply  this  method  for  a  source  signal  such  as  a 
gated  sinusoid.  More  importantly,  the  delay  estimates  ob¬ 
tained  via  this  method  are  biased  and  do  not  correspond  to 
the  true  parameter  estimates,  as  demonstrated  in  Sec.  IV  A. 
In  addition,  contiguous  frequency  samples  have  to  be  con¬ 
sidered,  which  might  not  be  desirable,  as  will  be  explained 
later. 

All  the  above  methods  require  the  source  signal  to  be 
known.  In  this  paper,  we  develop  an  algorithm  for  the  case  of 
a  known  source,  and  then  extend  this  algorithm  to  the  case  in 
which  the  source  signal  is  not  known  precisely,  but  is  as- 
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sumed  to  belong  to  a  parametric  class  of  signals.  That  is,  we 
assume  that  the  transmitted  signal  waveform  depends  on  the 
values  of  certain  constant  parameters  {i.e.,  the  parameters 
are  not  random  variables).  Our  objective,  then,  will  be  to 
obtain  good  estimates  for  the  delays  and  at  the  same  time 
extract  the  parameters  of  the  source  signal.  An  example  of  a 
parametric  class  of  signals  is  the  class  of  rectangular  pulses. 
A  signal  in  this  class  is  characterized  by  three  parameters; 
duration,  amplitude,  and  starting  point. 

In  the  next  section,  we  set  up  the  least-squares  problem 
to  be  solved.  In  Sec.  II,  we  develop  an  algorithm  assuming 
that  the  transmitted  signal  is  known.  Section  III  extends  this 
algorithm  to  the  case  when  the  signal  is  unknown,  but  can  be 
described  by  a  set  of  parameters.  In  Sec.  IV  we  give  the  re¬ 
sults  of  the  known  signal  and  the  unknown  signal  algorithms 
using  simulated  as  well  as  experimental  data.  The  transmit¬ 
ted  signal  is  a  gated  sinusoid  in  these  examples.  In  Sec.  V  we 
give  a  discussion  of  the  results.  Finally,  Sec.  VI  contains  the 
conclusions. 


I.  THE  LEAST-SQUARES  ESTIMATOR 

Assuming  the  source  signal  to  belong  to  a  parametric 
class  of  signals,  the  sampled  received  waveform  can  be  mod¬ 
eled  as 

M 

/•(«]=  ^  0</i<jV  —  1, 

A  =  1 

(2) 

where  0  is  the  vector  of  parameters  that  characterize  the 
source  signal.  If  the  noise  w[n]  is  white  Gaussian,  then  the 
least-squares  estimator  is  also  the  maximum  likelihood  esti¬ 
mator.'^  The  squared  error  function  is  given  by 

2  r[n]-  ^  a*s[«-r*;0]  ,  (3) 

n  ss  0  L  /c  =  ! 

and  the  objective  is  to  minimize  E  with  respect  to  all  of  its 
arguments. 

In  the  above  time-domain  formulation  of  the  problem, 
we  are  restricted  to  values  of  that  are  integer  multiples  of 
the  basic  sampling  interval.  If  a  more  accurate  estimate  is 
required,  then  one  has  to  resort  to  interpolation.  This  incon¬ 
venience  can  be  avoided  by  transforming  the  problem  to  the 
frequency  domain,  where  can  take  on  a  continuum  of 
values. 

The  data  record  in  ( I )  exists  in  the  time  interval  [0,7] . 
If  this  data  record  is  set  to  zero  for  r  <  0  and  t>T,  then  the 
summand  in  Eq.  (3)  will  be  zero  outside  the  range 
0<.n<N  —  1 .  Thus  we  can  extend  the  summation  range  from 
minus  infinity  to  plus  infinity  in  (3)  and  invoke  Parseval’s 
theorem  to  get 

£'(0,r*,aj 

=:_L(  _5(«;0)  y  'dry,  (4) 

2rr  J  „  1 

where /?(<y)  and5(<y;0)  are  the  discrete-time  Fourier  trans¬ 
forms  ofrjn]  and5[n;0],  respectively. 

The  problem,  now,  is  to  approximate  /i  ( ry )  by  a  weight¬ 


ed  sum  of  complex  expuncmials  For  compuici  implciiicnia- 
tion,  we  sample  in  ilie  frequency  doni.un  ai  miei  sals  of  Ai  r 
Hence, 

e  ■  '  -e  ■  ■  c  ■  , 

where 

■<i  -- 

The  integrals  are  then  replaced  by  sum>  and  the  squared 
error  function  is  given  approximately  by  (after  omitting  the 
scale  factor  l/2;r) 

£(0,>I,,Ui)  ==  y  I2?  [nl  -  S[n;01  V  . 

rt  <)  *  k  I 

(6) 

where  A{/i;01  ~  5(nAiy;0)  and  /?(nA<y)  The 

above  error  expression  will  be  used  to  develop  both  the 
known  signal  and  the  unknown  signal  algorithms.  For  the 
known  signal  ca.se.  the  parameter  vector  0  in  Eq.  ( 6 )  is  fixed 
and  is  not  included  in  the  minimization,  whereas  for  the 
unknown  signal  case  Eq.  (6)  is  minimized  w  ith  respect  to  all 
the  parameters. 

11.  THE  KNOWN  SIGNAL  ALGORITHM  FOR 
NARROWBAND  SIGNALS 

In  this  section  we  consider  the  case  when  the  source 
signal  is  known.  In  the  next  section,  we  show  how  the 
known-signal  algorithm  can  be  used  iteratively  in  the  ca.se 
when  the  signal  is  not  known. 

When  the  signal  is  narrow  band,  most  of  the  energy  of 
the  signal  is  concentrated  in  the  passband.  For  example,  if 
the  signal  were  a  gated  sinusoid,  most  of  its  energy  is  concen¬ 
trated  near  the  main  lobe  in  the  frequency  domain.  Recall 
that  the  spectrum  ot  a  gated  sinusoid  is  a  sine  function  cen¬ 
tered  at  the  frequency  of  the  sinusoid.  In  this  case,  we  do  not 
have  to  include  all  frequency  points  in  the  minimization  but 
consider  only  those  values  near  the  main  lobe.  Going  even 
further,  one  need  not  again  include  all  frequency  points  near 
the  main  lobe,  but  only  those  which  have  significant  signal 
energy.  That  is,  one  can  threshold  the  signal  spectrum  and 
consider  only  those  points  which  are  above  this  value.  This 
has  the  advantage  of  working  only  with  those  frequency 
points  having  good  signal-to-noise  ratios.  Also,  since  the  re¬ 
ceived  and  the  modeled  signal  are  real,  only  half  the  signal 
spectrum  needs  to  be  considered  because  of  its  conjugate 
symmetry.  The  expression  for  the  error  in  Eq.  (6)  is  thus 
rewritten  as, 

£(A.,a)  =  ^  2il[n]— 5[n]  ^  a^e'^""  ,  (7) 

nr,  t  k  ■  ■  \ 

where  A.  is  the  vector  of  scaled  delays  [seeEq.  (5)]  and  a  is  a 
vector  of  amplitudes.  (The  tildes  are  used  now  to  avoid  cum¬ 
bersome  notation  in  subsequent  developments. )  The  nota¬ 
tion  ne  J  means  that  the  summation  is  over  those  values  of  n 
which  belong  to  the  set.  1  ,  and  not  necessarily  over  contigu¬ 
ous  values.  The  set ,  )  consists  of  those  frequency  points  at 
which  the  magnitude  of  the  signal  spectrum  exceeds  a 
threshold.  A  rule-of-thumb  is  to  set  this  threshold  to  roughly 
one-twentieth  to  one-tenth  the  peak  magnitude  of  the  signal 
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spedtrum.  Since  matrix  notation  leads  to  succinct  expres¬ 
sions,  let  us  define  the  following: 

/i, 

r=(R{q,]  /?(<?,]  •••  R[q,})', 

a  =  (ti,  a,  •••  u.v,)''. 


S  =  diag(S[q,]  Slq^] 

■■■  ■S(?,v]). 

II 

< 

^rl 

V"”'  •  ■  • 

P(X)  =SA(X). 

(8) 

where  9*6.  )'",  k—  1,2 . N  denote  the  frequency  samples 

that  are  selected.  Hence  Eq.  (7)  can  be  recast  as 

£(A.,a)  =  ||f-P(X)al|^  (9) 

This  is  the  error  expression  considered  in  Ref.  9.  However,  in 
that  work  the  amplitudes  were  allowed  to  be  complex.  As 
mentioned  in  the  Introduction,  this  leads  to  a  biased  estimate 
of  the  delay  parameters  if  the  SNR  were  not  sufficiently  high 
( we  show  this  in  Sec.  IV  A) .  But  an  important  advantage  of 
allowing  complex-valued  amplitudes  is  that  the  correspond¬ 
ing  error  surface  is  reasonably  smooth,  making  it  easier  to 
find  Its  global  minimum.  On  the  other  hand,  while  the  global 
minimum  of  the  error  surface  defined  by  Eq.  (9)  with  the 
amplitudes  constrained  to  be  real  yields  the  true  delay  pa¬ 
rameters,  finding  this  minimum  is  extremely  difficult  as  the 
surface  is  characterized  by  numerous  closely  spaced  local 
minima.  We  now  outline  a  procedure  that  leads  us  to  this 
global  minimum,  starting  with  the  error  surface  correspond¬ 
ing  to  the  complex  amplitudes. 

The  key  step  that  helps  us  to  achieve  the  transition  from 
the  biased  error  surface  (due  to  the  complex  amplitudes)  to 
the  true  error  surface  (due  to  the  real  amplitudes)  is  to  re¬ 
write  Eq.  (9)  explicitly  in  terms  of  real  and  imaginary  parts 
of  each  of  the  terms  involved,  and  add  a  penalty  term  to  the 
imaginary  part  of  the  complex  amplitude.  That  is,  let 


£■„(>., a) 


+  a||a,| 


(10) 


where 


r,  =  Re{f},  r,  =  Im{f}, 

P, -Re{P(A.)},  P,  =  Im{P(A.)}, 
a,  =  Re{a},  a,  =  Im{a}. 

Rewriting  the  right-hand  side  of  Eq  ( 10)  a*"  a  sing'c 
norm,  we  get 


(11) 


where 


It  is  easy  to  see  that  for  a  =  0,  the  error  function  in  Eq.  (11) 
allows  the  amplitudes  to  be  complex  since  there  is  no  penalty 
attached  to  the  imaginary  part.  As  a  increases,  the  penalty 
on  the  imaginary  part  increases,  and  in  the  limit  as  a  —  a; , 
the  amplitudes  are  constrained  to  be  purely  real.  Note  that 
the  model  for  the  channel  assumes  the  amplitudes  to  be  real. 

With  the  above  modification  of  the  error  expression,  the 
procedure  to  minimize  £„  (A.a)  is  as  follows:  Initially,  a  is 
set  equal  to  zero.  Then  £,)  ( is  minimized  with  respect  to 
X  and  a.  The  minimum  yields  a  biased  estimate  of  X,  but  one 
which  is  reasonably  close  to  the  true  value.  To  estimate  the 
true  delay  values,  a  is  increased  gradually  and  the  minimum 
of  £„  (X,a)  is  found  for  each  a.  The  minimum  of  (X,a) 
yields  the  final  parameter  estimates.  This  method  of  using  a 
quadratic  penalty  function  to  impose  a  nonlinear  constraint 
is  known  to  converge  to  a  minimum  of  the  true  error  sur¬ 
face.'^ 

For  any  fixed  X  and  a,  the  best  a  which  minimizes 
£„(X,a)  is  given  by 

a=  (P"P)  - 'P"r.  (13) 

Substituting  this  value  of  a  in  Eq.  ( 1 1 ;  we  get 

£„(X)  =  ||[I-P(P"P)  ■■’P"]r|p 

=  ||P^r||^  (I'r) 

where  P^  =  I  —  P(P^P) '  'P",  and  now  the  error  is  a  func¬ 
tion  of  the  delay  parameter  vector  X  only.  In  the  next  subsec¬ 
tion,  we  present  a  Gauss-Newton  approach  to  find  the  mini¬ 
mum  of  £„  (X). 


A.  A  Gauss-Newton  approach  to  minimizing  £„  (X) 

A  common  algorithm  for  minimizing  an  objective  func¬ 
tion,  which  is  expressed  as  a  squared  norm,  is  the  Gauss- 
Newton  algorithm. This  algorithm  uses  a  first-order 
perturbation  expansion  to  convert  the  nonlinear  minimiza¬ 
tion  problem  to  a  linear  one.  A  sequence  of  linear  problems 
are  then  solved  until  the  solutions  converge.  In  this  section, 
we  derive  the  formulas  needed  to  implement  the  Gauss- 
Newton  algorithm  for  the  error  function  in  Eq.  (11). 

We  begin  by  considering  the  value  of  the  error  function 
at  an  increment  AX  away  from  the  nominal  value  of  X,  i.e., 

£„(X-f  AX)  =  ||[I-P(X  +  AX)(P"(X  + AX) 

XP(X  +  AX))  'P"(X  + AX)]r||v 

(15; 

To  simplify  the  above  equation,  we  consider  a  first-order 
perturbation  expansion  of  P( X  -f  AX) ; 

P(X-f- AX)-P4- AP,  :16) 

where  means  “equal  to  first  order,”  and  AP  is  a  matrix 
which  is  a  linear  function  of  AX.  Substituting  this  expansion 
of  P(X  -p  AX)  into  Eq.  (15)  and  retaining  only  first  order 
terms  (see  Ref.  15  for  details)  we  get 
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£„  (A  +  AA)  =  |!PV  -  P‘(AP)(P"P)  '  'P^'r 

-P(P"P)  '(AP)"P'r||-,  (17) 

where  the  superscript  H  denotes  complex-conjugate  trans¬ 
pose.  Using  the  definition  of  P  from  Eq.  (12)  it  can  be  veri¬ 
fied  that 


where 

Q  =  diag{9,  q2  ••• 

AA=(AA,  tiAj  AA^,)^ 

AA  =  diag{AA}. 

Since 


we  can  pull  diag{  AA;AA}  to  the  end  in  the  second  and  third 
terms  on  the  right-hand  side  of  Eq.  (17),  and  further  sim¬ 
plify  to  get 

£„(A  + AA)  =  ||x-BAAlp,  (19) 

where 

X  =  P'r, 

B'  =  P^S-diag{(P"P)  -‘P^r} 

+  P(P"P)  “'•diag{S’‘P'r}, 

B  =  B'(;,l:Af)  +  \  ):2M). 

( Note  that  we  have  used  standard  MATLAB  notation  in  the 
definition  of  B,  which  is  derived  from  B'  by  summing  the 
first  Af  columns  with  its  second  M  columns.)  FromEq.  (19) 
we  can  solve  for  the  best  A  A  (in  the  least-squares  sense)  as 

AA  =  B*x 

=  (B"B)  -  'B"x.  (20) 

Now,  A  is  replaced  by  A  AA  and  the  procedure  is  repeated 
till  iiAAII^  v^f.  In  many  cases,  it  rr.iglu  be  desirable  not  to 
make  large  changes  in  A  to  avoid  oscillations.  If  so,  the  least- 
squares  solution  in  Eq.  (20)  is  replaced  by  a  constrained 
least-squares  solution,'^  where  ||AA||2  is  restricted  to  not  to 
exceed  a  specified  value,  say  S. 

An  initial  value  for  the  parameter  vector  A  must  be  esti¬ 
mated  before  (20)  can  be  used  to  obtain  improved  param¬ 
eter  estimates.  A  standard  procedure  for  generating  initial 
estimates  is  the  coordinate  descent  algorithm  which  is  de¬ 


scribed  in  Refs.  17  and  18.  We  used  the  coordinate  descent 
algorithm  on  £„  ( A)  given  in  ( 14)  with  a  -  0  to  obtain  ini¬ 
tial  e.stimates  for  the  known  signal  algorithm. 


Hi.  THE  UNKNOWN  SIGNAL  ALGORITHM 

We  now  address  the  problem  of  time-delay  estimation 
when  the  source  signal  is  unknown.  In  this  case,  the  source 
parameter  vector  0  will  also  enter  the  miniini?.ation  in  Eq. 
(3).  In  the  algorithm  outlined  below,  we  assume  the  source 
to  belong  to  a  parametnc  class  of  signals,  viz.,  gated  sinu¬ 
soids  of  unknown  frequency  and  duration. 

First,  we  observe  that  the  source  signal  can  completely 
be  specified  (except  for  a  real-valued  scale  factor )  by  its  fre¬ 
quency,  duration,  starting  point,  and  phase.  If  the  frequency 
and  duration  are  such  that  there  are  many  cycles  of  the  signal 
present,  then  we  need  not  precisely  determine  the  phase, 
which  can  be  set  equal  to  zero.  Also,  since  all  time  delays  are 
relative  to  the  source  pulse,  the  starting  point  can  be  as¬ 
sumed  to  be  at  r  =  0,  without  loss  of  generality  . 

With  the  above  assumptions,  the  overall  problem  re¬ 
duces  to  minimizing 

min  y  i/?  [n] , 

Km. /.dot'll  err  I  I 

where  /and  d  are  the  frequency  and  duration  of  the  transmit¬ 
ted  signal,  respectively.  This  expression  is  not  easy  to  mini¬ 
mize  with  respect  to  all  the  parameters  simultaneously.  In¬ 
stead,  we  resort  to  an  iterative  minimization  technique 
outlined  below. 

(1)  Obtain  initial  estimates  of/^’  and  d'\  i.e.,  of  the 
unknown  frequency  and  duration,  respectively. 

(2)  Use/'  and  d '  in  the  known  signal  algorithm  to  esti¬ 
mate  A'  and  a'. 

(3)  Using  the  estimated  A'  and  a',  calculate /' '  '  and 
d'*\ 

(4)  Check  for  convergence  to  return  to  step  2. 

We  now  show  how  to  calculate /' "  '  and  d'  '  '  men¬ 
tioned  in  step  3.  Consider  the  error  function  for  a  given  A  and 
a: 

E{f,d)~  ^  R[n]  -  S[n-,f,d]  ^  . 

I  k  -  I 

Here,  E{  fd)  is  nonlinear  in/ and  d  and  an  analytical  mini¬ 
mization  is  again  very  difficult.  We  simplify  the  minimiza¬ 
tion  by  taking  advantage  of  the  fact  that  is  a  discrete  vari¬ 
able  since  the  signal  is  sampled,  i.e.,  only  finite  number  of 
values  fori/ are  possible.  Hence,  £(  f,d)  can  be  minimized  as 
follows: 

( 1 )  Carry  a  search  over  values  of  d  near  tne  previous 
estimate. 

(2)  For  each  d,  find  the  best  value  of / using  a  gradient 
descent  technique. 

(3)  Repeat  step  2  until  minimum  is  found. 

The  final  question  is  how  to  obtain /'  and  d‘\  In  our 
model  we  assume  that  the  first  path  is  the  least  attenuated 
one  and  that  its  amplitude  is  greater  than  that  of  the  noise. 
There  may  or  may  not  be  pulses  from  other  paths  overlap- 
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ping  with  the  first.  With  these  assumptions  we  get  rough 
initial  guesses  for/^  and 

(i)  Obtain  the  envelope  of  the  signal. 

(ii)  Obtain  the  high-amplitude  portion  of  the  received 

waveform  by  finding  that  part  of  the  signal  envelope  which 
exceeds  a  threshold.  Let  (t,  )  be  the  interval  in  which  the 

envelope  exceeds  the  threshold.  Choose  fi?“  =  (/,  —  )/2. 

( iii) /*’  is  obtained  using  any  standard  frequency  estima¬ 
tor  on  the  signal  in  the  interval  (r,  ./j ). 

This  completes  the  description  of  the  unknown  source 
algorithm. 

IV.  ALGORITHM  PERFORMANCE 

In  this  section  we  show  the  performance  of  the  algor¬ 
ithms  described  in  the  previous  sections  using  synthetic  sig¬ 
nals  as  well  as  experimental  data.  We  first  present  a  one- path 
example  to  illustrate  the  nature  of  the  and  (X) 

error  surfaces.  We  then  demonstrate  the  effect  of  increasing 
a  on  £„  (X).  Finally,  results  on  experimental  data  using  the 
known-  and  the  unknown-signal  algorithms  are  presented. 


A.  A  one-path  example  to  compare  error  surfaces 

As  mentioned  in  Sec.  II,  the  true  parameter  estimates 
are  obtained  from  the  global  minimum  of£,„  (X).  However, 
its  minimum  is  difficult  to  find.  On  the  other  hand,  the  error 
surface  £©  (X)  is  easier  to  work  with,  but  its  global  minimum 
yields  an  increasingly  biased  estimate  of  X  with  decreasing 
SNR.  We  now  illustrate  the  nature  of  £o(X)  and  £,^  (X) 
error  surfaces  with  a  one-path  example  to  highlight  the  diffi¬ 
culties  involved. 

A  synthetic  received  signal  is  considered.  It  was  con¬ 
structed  by  delaying  a  244-Hz,  40-ms  duration  sinusoid  by 
50  ms  and  adding  computer  generated  white  Gaussian  noise. 
The  error  surface  (plotted  as  a  function  of  relative  delay) 
should  have  a  minimum  at  /  =  0.05  s.  We  consider  two  cases, 
viz.,  high  and  low  SNRs.  (We  use  the  terms  "high”  and 
“low”  SNR  in  a  qualitative  sense. )  Figure  I  shows  the  high 
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FIG.  1.  Received  signal  at  high  SNR  for  a  synthetic  one-path  example. 
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FIG-  2.  Error  surfaCvS  £o  (^)  and  (A.)  for  l  he  high -SNR  received  signal 
of  Fig.  1. 


SNR  signal.  In  Fig.  2  we  show  the  corresponding  error  sur¬ 
faces,  £o(X)  and  £^  (X).  The  smooth  error  surface  is 
£o  (X),  while  the  sinusoidal  surface  is  £^  (X).  Both  have  the 
same  global  minimum,  but  £n(X)  is  easier  to  minimize, 
since  it  is  unimodal.  However,  when  the  SNR  is  not  high 
enough,  the  global  minima  of  these  two  surfaces  can  be  dis¬ 
tinctly  different.  Figure  3  shows  the  low  SNR  received  wave¬ 
form,  with  the  corresponding  error  surfaces  in  Fig.  4.  In  this 
example,  the  minimum  of  £o(X)  is  approximately  at 
t  =  0.049  s,  while  that  of  £^  (X)  is  still  att  =  0.05  s.  Itcanbe 
seen  that  if  the  time-delay  obtained  from  the  biased  error 
surface  is  used,  the  true  error  can  be  quite  large.  This  is  more 
so  in  the  M-path  case.  Hence,  one  has  to  minimize  £.,  (X)  to 
obtain  unbiased  parameter  estimates.  This  fact  was  first  not¬ 
ed  in  Ref.  19. 
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FIG.  3.  Received  signal  at  low  SNR  for  a  synthetic  one-path  exaniple 
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F{G.  4.  Error  surfaces  £„(k)  and  £,  (A.)  for  the  low-SNR  received  signal 
of  Fig.  3, 


B.  The  effect  of  increasing  a 

In  order  to  obtain  unbiased  estimates,  we  need  to  go 
from  the  minimum  of £0  (A.)  to  the  true  minimum,  i.e.,  that 
of  (A.).  This  is  achieved  by  incrc''.sing  a.  In  Fig.  5,  the 
smooth  error  surface  corresponds  to  a  =  0,  i.e.,  for  which 
the  amplitudes  are  complex.  As  a  increases,  the  norm  of  a, 
decreases  and  is  gradually  transformed  to  the  true 

error  surface,  as  can  be  seen  from  Fig.  5. 

C.  Performance  with  experimental  data 
f.  The  experiment 

Transient  data  were  gathered  in  the  Atlantic  Ocean  on  a 
bottom-mounted  receiver  in  780  m  of  water.  The  experimen¬ 
tal  geometry  is  shown  in  Fig.  6. 

The  acoustic  source  was  at  a  depth  of  40  m  and  transmit¬ 


FIG.  6  The  geometry  of  the  channel  used  for  the  experiments 


ted  a  244-Hz  gated  sinusoid  of  40-ms  duration.  The  source 
signature  was  recorded  using  a  hydrophone  mounted  on  the 
source.  The  signature  is  shown  in  Fig.  7.  The  pulse  was  trans¬ 
mitted  as  the  source  ship  drifted  over  the  bottom  receiver 
shown  in  Fig.  6.  The  horizontal  range  is  estimated  to  be  100 
m. 

The  ocean  bottom  is  characterized  by  a  thin  sediment 
layer  over  a  highly  reflecting  basalt  as  shown  in  Fig.  6.  The 
sediment  varies  in  thickness  from  0  to  20  m.  For  this  prob¬ 
lem,  a  10-m  sediment  thickness  was  chosen.  This  environ¬ 
ment  was  modeled  with  a  fast  field  program,  SAFARI.^"*' 
A  broadband  Gaussian  pulse  is  transmitted  in  the  model. 
The  model  predicts  four  paths  shown  in  Fig.  6  ( note  that 
Fig.  6  is  an  artist’s  sketch  of  the  four  paths).  Path  D  is  the 
direct  path,  path  DB  is  the  direct  path  reflection  off  the  ba¬ 
salt,  path  S  is  surface  reflected,  and  path  SB  is  the  path  re¬ 
flected  from  the  ocean  surface  and  the  basalt.  The  mode! 
ocean  impulse  response  is  shown  in  Fig.  8  with  the  four  paths 
labeled.  Pressure  release  surfaces  such  as  the  air/ocean  inter¬ 
face  cause  a  180°  phase  shift  in  the  reflected  signal,  causing 
the  negative  peaks  seen  in  Fig.  8.  The  received  signal  is 
shown  in  Fig.  9. 


xIO* 


Relative  Time  Delay  (secs.) 

FIG.  5.  Effect  of  increasing  a  on  the  constrained  error  surface  E„{K). 


Time  f  seen  J 

FIG.  7.  The  transmitted  (source)  signal:  a  244-H7.  gated  sinusoid. 
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FIG.  8.  Broadband  Gaussian  pulse  propagaied  with  SAFARI,  source 
depth  =  40  tn.  range  =  0. 1  km. 


2.  The  known  signal  algorithm 

Using  the  known  signal  algorithm  presented  in  Sec.  II, 
the  time  delays  and  amplitudes  were  estimated.  With  these 
parameters,  the  received  signal  was  reconstructed  (see  Fig. 
10).  The  residual  error  between  the  received  signal  and  the 
reconstructed  signal  is  shown  in  Fig.  1 1.  It  is  seen  that  the  fit 
is  very  good,  indicating  the  accuracy  o^  the  estimates.  The 
estimated  parameters  are  shown  in  Table  I  in  the  column 
labeled  “K.S.  estimate.”  Note  that  the  signs  of  the  estimated 
amplitudes  need  not  necessarily  agree  with  that  predicted  by 
the  channel  model  (see  Fig.  8),  which  shows  that  the  first 
two  amplitudes  are  positive  while  the  second  two  are  nega¬ 
tive.  The  signs  of  the  amplitudes  were  not  constrained  in  this 
algorithm;  however,  such  a  constraint  would  be  easy  to  im¬ 
pose  using  an  additional  penalty  function. 

3.  The  unknown  signal  algorithm 

Next  we  applied  the  unknown  signal  algorithm  outlined 
in  the  previous  section  and  estimated  the  parameters  of  the 


FIG.  9.  A  record  of  experimental  data  containing  a  received  signal  with  four 
overlapping  paths. 
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FIG.  10  The  reconstructed  received  signal  using  channel  parameters  esti¬ 
mated  by  the  known-signal  algorithm 


transmitted  signal,  as  well  as  the  delays  and  amplitudes  in 
the  received  signal.  These  estimates,  along  with  those  ob¬ 
tained  from  the  known  signal  algorithm,  are  presented  in 
Table  I. 

The  estimated  source  parameters  agree  quite  cio.sely 
with  the  nominal  values  (244  Hz  and  40  ms)  of  the  actual 
source  signal.  The  channel  parameters,  too,  agree  with  those 
obtained  by  the  known  signal  algorithm,  although  they  are 
not  in  the  same  order.  The  first  three  paths  identified  by  the 
known  signal  algorithm  correspond  to  paths  1,  3,  and  4 
found  by  the  unknown  signal  algorithm.  There  are  many 
different  parameter  sets  for  this  problem  which  all  give  rea¬ 
sonably  good  fits  to  the  data.  The  efficacy  of  the  unknown 
signal  algorithm  can  be  judged  by  reconstructing  a  signal 
using  the  estimated  parameters,  and  computing  the  residual 
error  between  the  received  and  the  reconstructed  signals. 
The  reconstructed  signal  is  formed  by  convolving  the  mod¬ 
eled  source  signal  (a  zero-phase,  245-Hz  sinusoid  of  43-ms 
duration)  with  the  estimated  channel.  Convolution  is  ac¬ 
complished  by  multiplication  in  the  frequency  domain  fol- 
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FIG.  1 1.  The  residual  error  (received  signal  minus  reconstructed  signal) 
for  the  known-signal  algorithm. 

Vaccaro  eta/.:  Least-squares  lime-delay  estimation 


216 


216 


TABLE  I.  Parameter  estimates  using  the  known  signal  and  the  unknown 
Sir  la!  algorithms. 


K.S.  e''’mate 

U.S.  estimate 

/(Hz) 

245 

/)  (ms) 

43 

r. 

0.1544‘ 

0.2547'’ 

~  -n, 

0.0360 

0.0198 

7)  -  r, 

0.0453 

0.0360 

U  -  T, 

0.0898 

0,0435 

Jy 

0.7913 

0.8287 

<^1 

0.4925 

-0.1597 

-0.1724 

0,4926 

0.0864 

0.1810 

‘  Relative  to  the  source  pulse  shown  in  Fig.  7. 

’’  Relative  to  estimated  source  pulse  starting  at  t  =  0. 


lowed  by  an  inverse  DFT.  The  corresponding  residual  error 
is  shown  in  Fig.  1 2.  The  fit  in  this  case  is  poorer  compared  to 
the  known  signal  algorithm.  The  reason  for  this  is  explained 
in  the  next  section. 

V.  DISCUSSION 

The  biased  error  surface  corresponding  to  the  complex 
amplitudes  is  easier  to  minimize  as  it  is  reasonably  smooth. 
For  the  one-path  example  of  Sec.  IV  A,  th  "surface  was  seen 
to  be  unimodal,  i.e.,  has  no  local  minima.  In  higher  dimen¬ 
sions,  the  surface  is  still  reasonably  smooth  as  opposed  to  the 
«  =  00  case,  but  is  no  longer  unimodal.  Hence,  one  has  to 
begin  any  minimizing  routine  with  a  reasonable  initial  esti¬ 
mate  for  k  to  reach  the  minimum  of  E^lk),  to  avoid  getting 
stuck  in  a  local  minimum.  We  have  found  that  the  coordi¬ 
nate  descent  algorithm  provides  good  initial  estimates.  The 
error  surface  being  smooth  for  a  =  0  is  not  true  for  an  arbi¬ 
trary  signal.  If,  for  example,  the  source  signal  consisted  of 


FIG.  12.  The  residual  error  {received  signal  minus  reconstructed  signal) 
for  the  unknown-signal  algorithm. 
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FIG.  13.  The  residual  error  between  the  source  signal  of  Fig.  7  and  the 
optimum  zero-phase  sinusoid  parametric  model. 


multiple  sinusoids,  even  for  a  =  0  the  error  surface  would 
possess  numerous  local  minima.  We  are  currently  develop¬ 
ing  a  technique  which  will  overcome  this  problem  by  parti¬ 
tioning  the  frequency  axis  into  bins.  An  algorithm  similar  to 
that  given  in  this  paper  is  used  in  each  bin. 

The  rather  large  residual  error  obtained  from  the  un¬ 
known  signal  algorithm  is  now  explained.  This  poor  fit  is  not 
due  to  any  inadequacy  in  the  unknown  signal  algorithm  per 
se,  but  due  to  the  fact  that  the  actual  source  could  not  accu¬ 
rately  be  modeled  as  a  perfect  sinusoid.  To  demonstrate  this, 
we  ran  the  unknown  signal  algorithm  by  considering  the 
recorded  source  to  be  the  received  signal  and  tried  to  model 
it  by  a  sinusoid  with  M  —\.  The  estimated  frequency  and 
duration  for  the  source  were  found  to  be  f—  243  Hz  and 
af  =  44  ms,  respectively.  The  residual  error  is  shown  in  Fig. 
13.  There  is  considerable  mismatch  at  the  beginning  and  at 
the  end  of  the  source  pulse  due  to  gradual  signal  build  up  and 
decay.  This,  therefore,  is  the  reason  for  the  spikes  in  the 
residual  error  of  Fig.  12,  wherein  the  unknown  signal  algo¬ 
rithm  assumed  a  perfectly  sinusoidal  source. 

VI.  CONCLUSIONS  AND  FURTHER  WORK 

The  residual  error  for  the  experimental  data  of  Sec.  IV  C 
using  the  known  signal  algorithm  was  quite  small.  We  tried 
the  same  procedure  on  a  different  experimental  data  set 
which  used  a  chirp  signal  as  the  source  pulse.  The  bottom 
reflecting  surface  was  very  rough  in  this  case.  The  known 
signal  algorithm  was  unable  to  improve  the  fitting  error  be¬ 
yond  a  certain  point.  There  still  appeared  to  be  considerable 
signal  structure  in  the  error  residual,  indicating  scope  for  a 
better  fit.  We  feel  that  this  might  be  due  to  inadequacy  of  the 
channel  model  given  in  Eq.  ( 1 ) ;  further  work  is  being  earned 
out  to  verify  this. 

To  summarize,  in  this  paper  we  pointed  out  the  impor¬ 
tance  of  constraining  the  amplitudes  to  be  real  But  the  error 
minimization  in  this  case  turned  out  to  be  difficult.  The  com¬ 
plex  amplitude  case  is  easier  to  solve  but  yields  biased  pa- 
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rameter  estimates.  We  then  presented  an  algorithm  which 
finds  the  minimum  of  the  true  error  surface  starling  with  the 
biased  error  surface.  When  initial  parameter  estimates  are 
obtained  using  coordinate  descent,  our  algorithm  finds  the 
global  minimum  of  the  error  surface  in  spite  of  the  existence 
of  numerous  local  minima.  The  source  was  assumed  to  be 
known,  initially.  Finally,  even  if  the  source  were  unknown,  it 
was  shown  that  it  can  be  estimated  if  it  belongs  to  a  paramet¬ 
ric  class  of  signals. 
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